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The study of fluctuations in gene regulatory networks is extended to the case of Gaussian colored 
ON ' noise. Firstly, the solution of the corresponding Langevin equation with colored noise is expressed in 

terms of an Ito integral. Then, two important lemmas concerning the variance of an Ito integral and 
the covariance of two Ito integrals are shown. Based on the lemmas, we give the general formulae 
for the variances and covariance of molecular concentrations for a regulatory network near a stable 
\ equilibrium explicitly. Two examples, the gene auto-regulatory network and the toggle switch, are 

presented in details. In general, it is found that the finite correlation time of noise reduces the 
^"*^ ■ fluctuations and enhances the correlation between the fluctuations of the molecular components. 

CN \ PACS numbers: 87.16.Yc, 87.10.-e, 05.40.a 

I. INTRODUCTION 

A regulatory network of gene expression consists of a group of genes which co- regulate one another's expression. 
The networks provide a fundamental description of cellular function that is the subject of the recently emerged field 
of systems biology [l[ . The advance of experimental techniques in constructing synthetic networks provided with the 
basic elements, such as a switch @, [H, H| and an oscillator [H, @, for the design of biological circuits For such 
elements, the characteristic features are determined mainly by the regulation scheme, and fluctuations always appear to 
be significant due to low copy numbers of molecules and stochastic nature of biochemical reactions 0, Q ■ In clarifying 
^v^j \ the nature of regulation mechanism, one of the important questions is to understand the way of accommodating 
the fluctuations meanwhile maintaining the stability for a genetic network. Recently, significant progress has been 
made along this aspect. One of the noticeable examples is the auto-regulatory network of a single gene for which, 
. the protein encoded in the gene serves as the regulator of itself through either negative or positive feedback. Such 
autoregulation is an ubiquitous motif in biochemical pathways yj. It was demonstrated by Becskei and Serrano that 
an autoregulatory network with negative feedback may gain stability Q. Further analysis given by Thattai and van 
• Oudenaarden [10( and by Ozbudak et al. indicate that noise is essentially determined at the translational level and 
negative feedback can suppress the intrinsic noise. Another example is the genetic toggle switch. Such a switch consists 
of two transcription factors which regulate each other's synthesis negatively [ij . Detailed analysis given by Cherry and 
Adler [l2| shows that the cooperative binding of two or more proteins in negative regulation is required, in general, for 
a switch to have two distinct stable states. Further study carried out by Warren and ten Wolde [13| indicates that the 
switch stability can be enhanced by overlapping the upperstream regulatory domains such that competing regulatory 
molecules mutually exclude each other, and robustness against biochemical noise may provide a selection pressure that 
drives operons together in the course of evolution. Recently, the results reported by Loinger et al [14[ reveal that 
a suitable combination of network structure and stochastic effects gives rise to bistability even without cooperative 
binding. 

Stochastic fluctuations associated with a system are often assumed to be Gaussian white noise in nature. However, the 
zero correlation time for white noise assumes an infinite relaxation time. Thus, it is important to incorporate the effect 
of finite correlation time of noise into the study of stochastic fluctuations. As the treatment for the effect induced by the 
Gaussian colored noise in a regulatory network being still lacking in the literature, we intend to fill up this gap in this 
paper. In modelling the dynamics of a regulatory network, rate-equation approach is often used; the approach reflects 
the macroscopic observation of deterministic nature. Noise-induced effects may be incorporated into the framework by 
employing the master equation and then proceeding via stochastic Monte Carlo simulations [l5j as was done for the 
first time for gene regulatory networks by Arkin et al 0|. In general, master equations are discrete in nature. By 
using the technique of 57-expansion [l7j , one can convert a master equation to continuous Fokker-Planck equation which 
then can be treated analytically by various approximations. Based on the Fokker-Planck equation, Tao and Tao et al. 
derived the corresponding linear noise Fokker-Planck equation which is suitable for the study of fluctuations caused by 
the white noise [3 [l^| . One may include the Gaussian colored noise into the frame of Fokker-Planck equation by using 
the scheme of unified colored-noise approximation [20| or via kinetic Monte-Carlo simulations as was implemented for 
gene regulatory networks in '21] ■ Alternatively, we first establish the equivalent Langevin description for a system with 
colored noise in the linearized region of a stable point [22| . The solution of the corresponding Langevin equation is given 
in the form of an Ito integral, and the fluctuations of molecular concentrations then can be evaluated by using the two 
lemmas, which concern with the variances and covariances of Ito integrals, shown in this work. As the correlation time 
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of noise is set to vanish, we recover the results of fluctuations for white noise. Based on this approach, we analyze the 
stochastic fluctuations of autoregulatory networks and toggle switches with general cooperative binding, and the results 
are found to be in good agreement with those obtained from numerical simulations. In general, the appearance of finite 
correlation time of noise decreases the fluctuations of a system and enhances the correlation between fluctuations. 

This paper is organized as follows. In section II, we first set up the model system for the study of the steady-state 
statistics of a regulatory network near a stable equilibrium. The rate equations are first linearized about the stable 
equilibrium. Starting with the rate equations, we then use the technique of ^-expansion to obtain the corresponding 
Fokker-Planck equation of the system. In section III, the equivalent Langevin description of the Fokker-Planck equation 
is given, and the solution of the Langevin equation with colored noise is expressed in terms of an Ito integral. Then, two 
important lemmas concerning the variance of an Ito integral and the covariance of two Ito integrals are shown. Based 
on the lemmas, we obtain the variances and covariance of molecular concentrations for a two-dimensional regulatory 
network near a stable equilibrium. In sections IV, the formulae resulting from the lemmas are applied to auto-regulatory 
networks and toggle switches to study the stochastic fluctuations of the systems. The comparison between the results 
and numerical simulations is given. In particular, the effects of noise correlation time on the amount of fluctuations and 
on the correlation between the two molecular components are analyzed. Finally, we summarize the results in section 
V. 



II. FOKKER-PLANCK EQUATION 



Consider a two-dimensional regulatory network of gene expression defined by the macroscopic rate equation, 



dx (t) = / (x) dt, (1) 

where the state variables x T (t) = (xi (t) ,X2 (t)) are molecular concentrations, and the forces f T (x) = (/i (x) , fa (x)) 
determine the time evolutions of the state variables. Here, the superscript r denotes the transpose of a vector. A 
stable point x* , specified by zero force / (x*) = 0, represents a stable stationary state. For the region near the stable 
equilibrium, the leading order of the drift force in x gives 

f(y)=F(x*)-y, (2) 

where the stable point is chosen as the origin, y — x — x*, in our two-dimensional space, and the elements of F (x*) 
are defined as (x*) = dfc (x) /dxj\ _ From now on, we drop the arguments whenever the matrix elements 
are understood as functions of the equilibrium stable point x*. The matrix F can be diagonalized by means of the 
transformation matrix P, 



F = P ■ 



-Ai 



P- 



-A 2 

where Ai, A2 > for the stable point, and det (P) = 1. The matrix P is set as 



P = 



P11 P12 
P21 P22 



and this gives the inverse as 



P- 



P22 
-P21 



~P12 
Pll 



In general, the drift force of Eq. (fT]) can be expressed as the sum of two terms, 

f(x)=R (x) - 6 • x, 



(3) 



(4) 



(5) 



(0) 



where the functions, R T (x) = (Pi (x) , R2 (x)), describe the synthesis or feedback regulation of molecule, and 6 is a 
2x2 constant matrix with the elements given by 8jj = 8ij6^ with the degradation rate 9^ 1 ' for molecular concentration 
Xj. Since the synthesis or feedback regulation of molecule concentration x, depends only on the concentration of other 
component, we have Ri (x) = Ri (xj) for j =/= i. The stochastic fluctuations can be incorporated into Eq. (fTJ) by 
means of the master equation approach. For this, we introduce the volume factor f2 to relate the molecular numbers 
n T = (m, 712) to the concentrations x T = (xi,X2) as n T = Qx T . In terms of molecular numbers n, the corresponding 
master equation for Eq. (fT]) with the drift force given by Eq. © is 



dP(n,t) 
dt 



i=l 



e®rii) P(n,t) 



2 

fi^Pi(x) [E t - 
»=i 



l]P(n,t), 



(7) 



where P (n, t) is the probability distribution, the step operators Ei± are defined as Ei±G (rii) = G {ni ± 1) for a function 
of molecular numbers G (n), and the fact, Ri (x) = Ri (xj) for j ^= i, is used for the term in the second sum. Then, the 
technique of fi-expansion [17| is employed to transfer the discrete process of Eq. (|7|) to a continuous process described 
by the Fokker-Planck equation, 

^ + V. J( M) = 0, (8) 

where (V) r = {d/dx\, 8/8x2), p{x,t) is the distribution density, and j (x,t) is the density current given as 

j (x, t) = f (x) p (x, t) ~ [D (x) ■ V] p (x, t) . (9) 

Based on Eq. ([7]), we obtain the diffusion matrix D (x) of Eq. © as Z)y (x) — Sijd^ (x), which takes the diagonal 
form with the diagonal elements given as 

d (i) (z) = ^ [R i (x) + 9^x i ] . (10) 

For the linear region specified by Eq. ([2|) we obtain the corresponding Fokker-Planck equation by expanding the 
density current J (x,t) of Eq. © around the stable point x* . The result reads 

^^ + V.Jz(y,*) = 0, (11) 

where Jl {y, t) contains only the leading order terms of J (x, t) in I/O, 

J L (y, t) = [F (x*) -y — D (x*) ■ V] p L (y, t) . (12) 



This leads to an Ornstein-Uhlenbeck process [17|, |22j, [23| in which, the drift force is linear and the diffusion is given by 
a constant matrix. 



III. EQUIVALENT LANGEVIN DESCRIPTION 

For the Fokker-Planck equation of Eq. (fTTj) , we have the equivalent Langevin description specified by the stochastic 
differential equation, 

dy(t)=F-y (t) dt + -q (t) dt. (13) 

From hereafter the stochastic fluctuations, described by the variables rj T (i) = (771 (t) ,772 (£)), will be assumed to be 
Gaussian colored noises. The two independent colored noises are specified by the differential equation, 

dr] (i) = -r -r)(t)dt + r ■ A • dW (t) , (14) 

where the constant matrices, T and A, and the vector dW (t) are defined as follows. The elements of the T matrix are 
given as = Sij (l/rW), where rW is the correlation time of the noise rji and the Kronecker delta 5ij is equal to 1 
for i = j and otherwise. The A matrix is related to the diffusion matrix evaluated at the stable point, D = D (x*), 
by A • A = 2D. Based on the form of the diffusion matrix of Eq. (fT0|) . we have Ay = Sijy/2d^ (x*) for the elements of 
A. Furthermore, the variables dW T (t) = (dw\ (t) , dw2 (t)) describe two independent Wiener processes. By rewriting 
dW (t) — > £ (t) dt, the conditions, (£ (t)) = and (£ (t) ^ T (s)) = S (t — s) with the Dirac delta function S (t), specify the 
Gaussian white noise. 

The solution of Eq. IT4]) can be expressed in terms of Ito integral as 

T] (t) = 7] (0) exp (-tr) + / [exp (- (t - s) T)} ■ T ■ A ■ dW (s) . (15) 
Jo 

This yields the correlation function of the Gaussian colored noise as 

(r)(t)r) T (t')} = T-Dexp(-\t-t'\T). (16) 
Then the solution of Eq. (Tl"3"l) can be written as 



y(t)=y (0) exp (tF) + / [exp (t - u) F] ■ \ rj (0) exp (-uT) + / [exp (- (u - s) T)} ■ T ■ A • dW (s) \ du. (17) 



Since we are interested in the fluctuations around the stable point x* , only the asymptotic behavior of the solution 
matters, and we may set the initial condition as y (0) = and 77 (0) = to obtain 



y(t) = J [exp (t -u)F]- |^ [exp (- (u - s) T)} ■ T ■ A • dW (s) } du. 
Moreover, the order of the double integrations can be changed properly to yield 



y(t) 



[J [exp (t —u)F]- [exp (- (u — s)T)]-T ■ Adw j • dW (s) . 



(18) 



(19) 



Thus, the form of the solution, in general, is an Ito integral. 

Two important lemmas for the evaluations of variances and covariances of Ito integrals are shown as follows. Consider 
an Ito integral in the form of 



1(f) 



V(t,s)dw(s) 



for a random Wiener process dw. The integral can be expressed as the discrete form, 

N-l 

I(t)= Urn ^2v(t,s k )[w(s k+1 )-w{s k )}, 



k=0 



with sq — and st — k (t/N). This leads to the variance of I (t) as 



N-l 



cr 2 [I{t)]= lim Y\ V 2 (t, s k ) a 2 [w (s k+1 ) - w (s k 

TV— >OG * — 4 



(20) 



(21) 



(22) 



k=0 



Here the variance of / (t) is defined as a 2 [I (t)] = (I 2 (t)) — (/ (t)) 2 with the expectation value (•) taken with respect 
to the distribution of noise at time t. Based on the equality for a Wiener process, tr 2 [w (s k +i) — w (s k )] = s k +i — s k , 
we can express Eq. f|22[) as the integral form and obtain the first lemma, namely, the variance of / (t) is 



a 2 [I(t)] 



V 2 (t,s) ds. 



(23) 



The first lemma can be further extended to the case of two dimensions in a straightforward way. Consider the Ito 
integrals, J T (t) = (J 1 (t) , J 2 (x)), defined as 



J (t) — / M (t, s) ■ dW (s) , 



where M (t, s) is a 2 x 2 matrix, 



M (t, s) 



mu(t,s) m 12 (t,s) 
m 2 i (t, s) m 22 (t, s) 



(24) 



(25) 



and the two variables dw\ and dw-i in dW T (t) = (dwi (t) , c£u> 2 (t)) describe two independent random Wiener processes. 
Following the first lemma along with the property (dW (t) dW T (s)) = 5 (t — s) dt, one can show that the variances of 
Ji (t) and J 2 (t) are 



^ 2 [Ji (*)] = / [ma (*, *) + m 2 i2 (t, a)] ds 



for i = 1,2, and the covariance is 



E [ Ji (t) J 2 (t)] = \ [m n (t,s)m2i(t,s) + mi2(t,s)m 2 2(t,s)] 
Jo 



ds. 



(26) 



(27) 



Here the covariance of J\ (t) and J 2 (t) is defined as E [J\ (t) J 2 (t)] = (Ji (t) J 2 (t)) — (Ji (t)) (J 2 (t)). These constitute 
the second lemma. 

The second lemma, Eqs. (|26p and (1271) . with the limit t — > 00 can be employed directly to determine the variances 
and covariance of y (t) of Eq. (I19|) near the stable equilibrium y = 0. The results thus obtained are summarized in the 
following. The variances of y\ and y 2 , referred as o\ and o 2 respectively, are 



a 2 (r(V (2) ) =d (1) 



4 



2g3ff4Cl 

Aicn A 2 c 2 i Ai + A 2 



dW 



2g 2 ? 2 

AiCi 2 A 2 c 22 Ai + A 2 



(28) 



and 



a 2 2 (VV (2) ) =d (1) 



4 



4 



^ 2gjci 

Aicn A 2 c 2 i Ai + A 2 

and the covariance between y\ and y^ denoted as E, is 

9293 . 9294 (9293 + <7294) Ci 



(rW T (2) ) =d W 



Aicn A2C21 



Ai + A2 



4 



4 



29394C2 

A 2 c 2 2 A1C12 Ai + A 2 



9i94 



9i93 



(9l93 + 9i94) c 2 



A1C12 A2C22 



Ai + A2 



(29) 



(30) 



Here, the quantities qi are defined as q\ — P11P12, 92 = P21P22, 93 = P11P22, an d 94 = P12P21 with pij specified by the 
transformation matrix P of Eq. ([5]); and the quantities Ci are q = (1/ci,) + (l/c 2 i) with cy defined as cy = AjrW + 1; 
for i, j = 1 and 2. 

The stochastic differential equations specified by Eqs. (fTB")) and ([H]) reduce to 



dy(t)=F-y (t) dt + A-dW (t) . 



(31) 



for the case of white noise. Similarly to the case of colored noise, the solution of Eq. (|3"Tj) can be put in the form of Ito 
integral, and we obtain the variances and covariance of y\ and j/2 as 



af )2 = d« 



'4 



4 



49394 

Ai A2 Ai + A 2 



S 2 
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Ai A2 Ai + A2 



(32) 



J0)2 



'4 



4 



Ai + A 2 



d^ 



'4 



4 



Ai Ai + A2 



49394 



(33) 



and 



E (0) = d (D 



9293 . 9294 2 (q 2 q 3 + q 2 qi) 



Ai 
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Ai + A2 



d^ 



9i94 , 9i93 2 {q x q z + q^) 



A1C12 A2C22 



Ai + A2 



(34) 



where the superscript (0) in a and E is used to denote the result of white noise. These are exactly the same as those 
given by Eqs. (^HI) - (EOJ) but with = = 0. Following this, by setting either = or = in Eqs. ([2^|) - ([50|) 
we have the results for the case of Gaussian colored noise in one component and white noise in the other. 



IV. RESULTS AND DISCUSSIONS 



The above results are applied to analyze the fluctuations of two systems, including auto-regulatory networks and 
toggle switches. The fluctuations of the system near a stable point are analyzed by measuring the variance of xi in 
terms of the Fano factors z/,-, defined as 

and the covariance of xi and X2 in terms of the correlation coefficient i?i 2 , defined as 

r„ ( T W r (2) N ) = ^ T ' T > (36) 



Here at (t^\t^) = \J a 2 (tW,t( 2 )) is the standard deviation of the concentration Xi near a stable equilibrium. Note 

that we add the superscript, (0), to a quantity to refer to the results of white noises, i.e. = and = 0. 

The Fano factor is equal to one, i/j = 1, for a Poisson process. Based on this, we refer a process with Fano factor 
smaller than one as sub-Poissonian and a process with Fano factor larger than one as super-Poissonian. Thus, the 
sub-Poissonian deviates from the inherent randomness of Poisson process in a opposite way to the super-Poissonian, 
the former suppresses the occurrence probability of the large deviations from the mean value x* meanwhile the latter 
increases it. 



-Ik 



-(1) T (2)> 



(35) 




FIG. 1: The fluctuations for autoregulatory networks: For the case of white noise, the Fano factors, v\ (hollow circles, left 
vertical) and v^ 1 (black dots, right vertical), versus the concentrations of mRNA (lower horizontal) and protein (upper horizontal) 
at the stable points for different Hill coefficients j3 are shown in (a). For the case of Gaussian colored noise, the Fano factors, 
v\ (r' 1 - 1 , r*- 2 ') and (r^ ,T^ 2 'J , and the correlation coefficients of noises R12 fr^, t^ 2 'J versus the correlation times of noises, 

t (lower horizontal) and r' 2 ' (upper horizontal), at the stable points for different Hill coefficients /3 are shown in (b), (c), and 
(d), respectively. 



A. Auto-regulatory network 

The two variables of auto-regulatory network, x\ and X2, are referred to the concentrations of mRNA and protein, 
respectively. For the functions Ri (x) in the drift force of Eq. ([5]), we adopt the most common noise-attenuating 
regulatory mechanism, called negative feedback and described by Hill function, 

*i(sa) = ?, (37) 



to regulate the production of mRNA, and set 



1 + (x 2 /k d y 



R 2 (x 1 ) = k 2 x 1 . (38) 



Here fc m ax is the maximum transcription rate of mRNA, kd is the binding constant specifying the threshold protein 
concentration at which the transcription rate is half its maximum value, (3 is the Hill coefficient, and k% is the translation 
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FIG. 2: The Fano factors of toggle switch, (hollow circles) and (black dots), versus the Hill coefficient /3 are shown in 
(a) for the stable point A and in (b) for the stable point B. Note that while v± and v^ 1 are shown with the same vertical scale 
in (a), the left vertical scale is for vf^ and the right vertical scale is for in (b). 



rate of protein k 2 . Then, a stable equilibrium, x* , can be characterized by two conditions: 9^6^ — r\ (x 2 ) k 2 > and 
ff(X) -\-ffi 2 ) > with n (x 2 ) = dR\ [X2) /® X 2\ X2 = X *- Subsequently, one can use Bendixson's criterion to further conclude 
that there are no cycles and only one equilibrium exists [24| . 

We mainly follow Refs. [13, Gil] to specify the values of the parameters as follows. The half-lifetimes of mRNA 
molecules and proteins are set as 2 minutes and 1 hour, respectively; this leads to 9^ = (In 2) /2 and 9^ = (In 2) /60 
in the unit of (min) -1 . The average size of a burst of proteins, b = k 2 /9 ( - 1 \ is set as 10. This leads to k 2 = 5 (In 2). By 
using the fact that the protein concentration is about 1200 when (3 = (no feedback), we set fc max = 3 [25| . To study 
the effect of the strength of negative feedback on the fluctuation of the system, we vary the parameters (3 from 2 to 10, 
while the kd value is fixed as 800. 

We first consider the case of white noise. Based on Eqs. ([32|) -(p4 | . we obtain the z^and values for (3 ranging from 
2 to 10, and the results are shown in Fig. 1(a). The characteristic features revealed from the results are summarized 

as follows. The v% value is much larger than the corresponding v[ ^ value, and it decreases as f3 increases. This leads 
to the well-known conclusions that the stochastic fluctuations occur mainly at the translation level, and the negative 
feedback may enhance the stability of the system. On the other hand, the value decreases from vf^ = 0.9844 at 

(3 = 2, reaches the minimum = 0.9827 at f3 = 4, and then increases to = 0.9907 at [3 = 10. All the i/[ ^ values 
are very close to but always less than one. This implies that the process of transcription is sub-Poissonian. We then 
consider the effect of non-zero correlation time of noise on the fluctuation of the system based on Eqs. (|2"8]) - (f3"0"|) . First, 
we show the plots of v x (t^,t^) and v 2 (t^,t^) versus (3 for (t^\t^) = (0.1,6), (0.2,12), (0.4,24), and (0.6,36) 
in Figs. 1(b) and 1(c). Here the values of (t^,-^ 2 )) are set as the multiples 10, 20, 40, and 60 of the half-lifetimes 
of mRNA and protein, respectively. The results indicate that the fluctuations are reduced by the amount proportional 
to the correlation time value of the noise. However, the reduction is less significant for protein owing to its long half- 
lifetime when compared with that of mRNA. Moreover, the effect of reduction, in general, is enhanced for a larger Hill 
coefficient, but this is not very significant for mRNA owing to the nature of the sub-Poissonian process. We also show 
the effect of the correlation time of noise on the correlation coefficient between x\ and x 2 in Fig. 1(d), where the plots 
of i?i2 (t^'jtC 2 ') versus (3 for five sets of (r^ 1 ',^ 2 ') values are displayed. The results indicate that the correlation 
coefficient R\ 2 increases with the correlation times of the noises but decreases with the Hill coefficient. Moreover, as 
the consequence of sub-Poissonian distribution for the mRNA component, the R\ 2 values between two components are 
very small, ranged between 0.0202 and 0.1213. 

The analytical results given by Eqs. (|2"5|) - (|3"0")) are obtained with the approximation of linearization about a stable 
equlibrium. To test the validity of such an approximation, we solve the stochastic differential equations numerically 
by using Heun's method. This numerical method is a stochastic version of the Euler method, which reduces to the 
second-order Runge-Kutta method in the absence of noise [26l | . The numerical simulations are in a very good agreement 
with the analytical results as shown in Figs. l(b)-l(d) for (3 = 2.0. 
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FIG. 3: The fluctuations of toggle switch for the case of Gaussian colored noise: The Fano factors, v\ (r^,T^) and ^2 (r' 1 ' , t' 2 ') , 
versus the correlation times of noises, r' 1 ' and r' 2 ' for different Hill coefficients (3 are shown in (a) and (6) for the stable point A 
and (c) and (d) for the stable point B. Note that r^ 1 - 1 is set to be equal to r^ 2 ' in the calculations of v\ and and the values 
of i^i and ^2 for r' 1 ' = r' 2 ' = are the results for the case of white noise. 



B. Toggle switch 



A toggle switch consists of two transcription factors with concentrations x\ and X2; the transcription factors can 
regulate each other's synthesis through the negative feedback mechanism. Such genetic circuits often exhibit two 
stable states, referred as A(x*,X2) and B (x\, for which, the x\ component is dominant with almost vanishing 
X2 component for the state A, and vice versa for the state B. The distinct two stable states can be switched either 
spontaneously or by a driven signal. A plasmid of this type in Escherichia coli has been engineered by Gardner, Cantor, 
and Collins Q. For this, correspondingly, the degradation rates are rescaled to 9^ = 8^ = 1; the regulation functions 
in the drift force of Eq. © then become 



Ri (a*) 



at 



.0 



(39) 



and 



R2 (an) 



«2 



1 



(40) 



Here the parameter values are ct\ — 156.25, ct2 — 15.6, j3 = 2.5, and 7 = 1.0. Such a specification leads to two stable 
states, A (155.7634, 0.0995) and 5(0.3324,11.7080). 

The system with the regulation functions of Eqs. (|39p and ([4^)1 exhibits the bistability over a wide rage of parameter 
values. In this study, we intend to analyze how the characteristics of the fluctuation change with the cooperative 
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FIG. 4: The correlation coefficients of noises of toggle switch, R12 (t' 1 ', t* 2 '), versus the correlation times of noises, r' 1 ' and 
for different Hill coefficients /3 are shown in (a) for the stable point A and in (6) for the stable point B. Here r' 1 ' is set to be 
equal to t^ 2 \ and the results for r' 1 ' = r' 2 ' = correspond to the case of white noise. 



binding of the system. Thus, we calculate the variances and covariances along the parametric path increasing the Hill 
coefficient j3 from 2.5 to 7.5 and keeping the other parameters to be the same as the previous values. As the (3 value 
varied from 2.5 to 7.5, the loci of the two stable states are changed as follows. The x\ value increases slightly from 
155.7634 up to 156.2500 and the x* 2 value decreases insignificantly from 0.0995 down to 0.0992 for the stable state A; 
meanwhile, the x\ value decreases from 0.3324 down to 1.759 • 10~ 7 and the x* 2 value increases from 11.7080 up to 
15.6000 for the stable state B. Accordingly, increasing the (3 value will enhance the major component and suppress the 
minor component, and it will shift the location of the stable state B more significantly than that of A. 

We first consider the case of white noise for the study of fluctuations. The results of i^and v% for different (3 values 
are shown in Figs. 2(a) and 2(b) for the stable states A and B, respectively. As indicated by the numerical results, the 
2/^ and ^2°^ values all are larger than one; thus, the fluctuations are caused by super-Poissonian processes for both A 
and B. This is opposite to the case of the auto-regulatory network, and it agrees with the results obtained by Tao [Lsl ]. 
Moreover, the system in the stable state B always possesses a larger deviation from Poissonian than that for A. But, the 
maximum deviation occurs at [3 = 2.5 for both A and B at which, we have = 1.0517 and — 1.0042 for state A 

and = 2.0580 and = 4.7268 for state B. Note that there is a big drop in the j^and v^p values for the system 
in the state B when the (3 value increases from 2.5 to 3.0. As the values of the finite correlation time of the noises set in, 
the resultant values of v\ (t^,-^ 2 )) and i>2 (t^,-^ 2 )) for the values of r' 1 ' = t^ 2 ) ranging from to 1 with (3 = 2.5, 3, 
and 6 are shown in Figs. 3(a) and 3(b) for the states A and in Figs. 3(c) and 3(d) for the state B. The results indicate 
a similar feature as that for the case of the auto-regulatory network, namely the fluctuations are reduced by the noises 
being correlated, and the longer is the correlation time, the bigger amount the fluctuation decreases. We also show the 
plots of R12 (t^,-^ 2 )) versus t^ 1 ' (= t^ 2 )) for three different (3 values in Figs. 4(a) and 4(b) for the stable states A 
and B, respectively. All R12 (tW,t( 2 )) values shown in the figures are negative. This implies that the fluctuations of 
two components are anti-correlated, which reflects the fact that the two components are negatively regulated with each 
other. In particular, the system in the state B with [3 = 2.5 and 3.0 is highly anti-correlated with = —0.8070 and 
—0.4453, respectively. Also, the anti-correlation of fluctuations is enhanced as the finite correlation times of noises are 
set in. 

The results shown in the above are also compared with those obtained from numerical simulations. The comparison 
is shown in the insets of Figs. 3(a)-3(d) and Figs. 4(a)-4(b) for {3 — 2.5, and much larger differences in the Fano factors 
have been observed for the stable point B than that for the stable point A. 



V. SUMMARY 



Noise due to stochastic fluctuations is always present and essential in the gene expression process due to low copy 
numbers of molecules and stochastic nature of biochemical reactions. A framework has been set up for the study of the 
noise by constructing the equivalent Langevin description of the system. By taking then the fluctuation as Gaussian 
colored noise, we first solve the Langevin equation, and then obtain the general formulae for the variances and covariance 
of a two-dimensional system near a stable equilibrium based on Ito calculus. For the latter, two important lemmas, 
concerning the variance of a Ito integral and the covariance between two Ito integrals, are established. 



We apply the general formulae to the auto-regulation network and the toggle switch for the study of stochastic 
fluctuations of the systems, in particular, the effect caused by finite correlation time. In general, as correlation time 
of noise is set in, it will decrease the fluctuation and enhance the correlation of fluctuations. For the auto-regulation 
network, the fluctuations of mRNA concentration mirror a sub-Poissonian process with the Fano factor very close to 
but less than one. Consequently, this leads to the correlation coefficient between the fluctuations of mRNA and protein 
concentrations to be very small. For the toggle switch it is peculiar for one of the two stable states when the system 
has the Hill coefficient (3 being around 2.5; the fluctuations are caused by a process far away from Poissonian, and the 
fluctuations between two molecular components are highly anti-correlated. Aside from this, the fluctuations are caused 
by super-Poissonian processes with Fano factors very close to but larger than one for both concentrations in any of 
the two stable states; and this then leads to small but negative correlation coefficients between the fluctuations of two 
components. Moreover, as it has been shown, these analytical results are in a very good agreement with the numerical 
simulations. 

To summarize, let us stress that it is more realistic to take into account the noise finite correlation time in order to 
study of stochastic fluctuations of gene regulatory networks. Meanwhile, in this work the noises of two components are 
still assumed to be independent, although the components couple with each other in the rate equations, ft was shown 
[27( 1 that increasing the coupling strength between two noises may drive the system to transit from a bistable stationary 
probability distribution to a mono-stable one. Thus, it would be interesting to study the effect of coupled noises on the 
stochastic fluctuations in the gene regula tory networks. Besides, it will be also very interesting to extend the work to 
the case of a non-Gaussian colored noise 128]. 
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